A non-Markovian model of rill erosion 

Michael Damron* CL. Winter^ 
December 2008 



Abstract 

We introduce a new model for rill erosion. We start with a network 
similar to that in the Dynamical Discrete Web** and instantiate a dy- 
namics which makes the process highly non-Markovian. The behavior 
of nodes in the streams is similar to the behavior of Polya urns with 
time-dependent input. In this paper we use a combination of rigorous 
arguments and simulation results to show that the model exhibits many 
properties of rill erosion; in particular, nodes which are deeper in the 
network tend to switch less quickly. 

1 Introduction 

1.1 Reinforcement and Rill Erosion 

Stochastic processes with reinforcement are inherently non-Markovian and there- 
fore may model some real phenomena more accurately than can their Markovian 
counterparts. Reinforcement is a mechanism that provides a bias to a system, 
making it more likely to occupy states the more often those states are visited. 
Some well-studied examples include variations on the urn of Polya (the original 
introduced in jl] and this and subsequent models studied, for example, in [1 
and [3]) and reinforced random walks ([3J US])- The infinite memory exhibited 
in these examples can force a system to spend most (or almost all) of its time 
in a small subset of its state space. Many natural phenomena exhibit similar 
behavior; for instance, the overall pattern of erosion on a hillslope is relatively 
stable once it is established, although small details of the pattern may change 
frequently and catastrophes that permanently alter it may occasionally occur. 
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We investigate a discrete time, infinite-memory random process defined on 
the nodes and edges of an oriented diagonal lattice (Figure [TJ that we propose 
as a simple model of hillslope erosion. The lattice starts out smooth in the sense 
that it has no edges initially, but it sprouts edges everywhere the instant the 
process starts, much as rain can start soil erosion everywhere on a hillslope at 
once. Edges may connect an interior node to two, one, or neither of the two 
nodes directly above it. Exactly one edge descends from each interior node, and 
it points either left or right. At every node and at every time step a simple two 
parameter reinforcing law, based on the entire history of the network above a 
given interior node, randomly determines the direction of the nodes descending 
edge and then is updated. Obvious modifications of these statements apply to 
nodes at the top or bottom (if one exists) of the lattice. 

The current pattern of connections among nodes represents the present state 
of the process, and the patterns stability measured by the tendency of the same 
state, or one similar to it, to occur on subsequent iterations of the process - 
represents the patterns strength as a memory. The degree of reinforcement is set 
by tuning two parameters, r and a. At any given moment the current pattern 
is a collection of dendritic networks that appears similar to drainage networks 
found in nature; indeed, lattice models have often been used to investigate the 
morphology of natural drainage networks (e.g. [T7]). We focus on the surficial 
dynamics of rill networks |10j . rather than their morphology. Put in terms of 
erosion, we are more interested in the process of erosion than we are in the 
result. 

The analogy between our model and erosion, specifically rill erosion, is 
straightforward: r can be interpreted as a rainfall rate (or cquivalcntly, as the 
rate of sediment generation) and a -1 as the resistance of soil to erosion, while 
the reinforcement dynamics correspond to the overland flow of water and sedi- 
ment down a hill. Rills are small, ephemeral channels that transport sediment 
down hillslopes when it rains |18] . They form when rainfall and runoff dislodge 
particles from the soil surface and transport them along flow paths governed 
by variations in the surface roughness of soils and the soil's ability to resist 
erosion. Flow depths in rills are typically on the order of a few centimeters or 
less, while the longest channels in rill networks can be several meters long. Pro- 
cesses affecting rill erosion take place over timcscalcs ranging from milliseconds 
to hours. 

The topology of rill networks is relatively unstable when compared to larger 
scale natural drainage systems (of which rills may be a part) like gulley systems 
and river basins. Rill networks are most unstable at their tops where bound- 
aries between rills and inter-rill areas are not well defined and shift often, but 
connectivity can change downhill as well, usually at a slower rate than uphill. 
Some rills grow throughout a rainfall event, others are filled by sediment and 
disappear, still others alternate. A detailed description of rill erosion 1) must 
account for complicated interactions among rainfall, soil properties, and topog- 
raphy 19J, and 2) often depends on obtaining a set of physical parameters that 
are difficult to measure. 

Despite the high degree of complexity of rill erosion at small scales, at macro- 



2 



scopic scales it is principally determined by particle detachment, overland flow, 
and sediment transport ( |16j). In turn, each of flow, detachment, and transport 
depends critically on the rate of rainfall and the soils resistance to erosion. It is 
not completely surprising that our simple two parameter model exhibits some 
important elements of the macroscopic behavior of rills formation. In fact, sim- 
ilar to rill erosion, each node in the model network switches direction infinitely 
many times but the switching rate depends on position up or down hill. Fur- 
thermore, floods that carry unusually large amounts of water and catastrophes 
that significantly alter the flow pattern occur occasionally in the model, as they 
do in nature. 

1.2 Definition of the Model 

Consider the vertices in the even sub-lattice of Z 2 which have second coordinate 
non-positive. That is, the set 

Z e „ en = {(#, y) £ 1? : x + y is even and y < 0} 

and edges 

Kven = {< ( x i V)> ( x + !> V - !) > : ( x > V) e Z Len} 

U{< (x,y),(x-l,y-l) >: (x,y) £Z 2 even } 

Let v = (x,y) be a node with left parent w\ = (x — l,y + 1), right parent 
w 2 = (x + l,y + 1) (for those nodes with second coordinate 0, parents will not 
exist), and with left child (x — l,y — 1), right child (x + l,y — 1). We will 
use the term depth k to refer to those nodes with y coordinate equal to 1 — k. 
Conversely, for any node v the term depth(w) will denote the numerical value 
of the depth of v. 

First we describe the algorithm for the behavior of v heuristically. At the 
end of the th second, v receives I v (0) = r units of rain (but does nothing else). 
During the n th second, for n > 1, the following sequence occurs: 

1. v flips a coin, heads-biased with probability P^(n), which reflects T^(n), 
the total "sediment" load v has sent to its left child by time n. 

2. If this coin shows heads (tails), v sends its current input of sediment 
I v (n— 1) to its left (right) child, v adds this number to the total sediment, 
T^(n) (T^(n)), it has sent to the left (right) for all time. 

3. v receives sediment load from 0,1, or 2 parents and receives r units of rain. 
Call the sum of these two I v (n). Increment time and return to step 1. 

The evolution of the node's behavior depends on two parameters: the rainfall 
rate r > and a term or 1 > that resists change. To make this precise, we 
make several definitions. We start by initializing variables. For each v € Zg„ en , 
let 
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2£(0)=r*(0)=T o (0)=0 

and 



l v (0)=r ,P£(0) = P?(0) = l/2 

For each n > 1, v € Zg„ en , we define a Bernoulli variable, D^(n), the biased 
coin, that is conditionally independent from vertex to vertex given the variables 
{D^(i) : v £ l? evenl i < n}, with parameter P^(n — 1). Next, let 

T v L (n) = I v (n - l)D^(n) + T v L (n - 1) 



T v (n) = I v (n - 1) + T v (n - 1) , T*(n) = T v (n) - T v L (n) 
We create the bias for the next coin: 

*<»> " ^TT? = (1 > 

T v (n) + 2a hjm + 2 V 

where i] = a/r compares the effect of the rain to the system's inherent resistance 
to change. In this paper, we shall always take r = 1 so that r\ — a. Last we 
define the input 



I v (n) 
and the filtration 



r depth(w) = 1 

I W1 (n - 1)(1 - (n)) + I W2 (n - l)D^ (n) + r otherwise 



T n - a({D^(k) : v G Z 2 even , k=l, n}) 

See Figure [T] for an illustration of the process at node v. 

Denote by d v = (d v (l), d v (2), ...) the sequence of directions that node v 
chooses (for example (L,R,L, ...)). At the end of time t, after all nodes have 
sent loads to their children, we may update certain edge variables. Define a 
sequence of edge configurations {t-u n } n >Q, where for each n, uj n is a map from 
Eg„ en — > {0, 1}, using the following rule. If the node v — (x,y) has d v (n) = L 
then let 



u n (< x- l,y- 1 >) = 1 , u> n (< x + l,y-l>)=0 
If, on the other hand, d v (n) — R then let 

u n (< x-l,y-l>) = , u n {< x + l,y-l>) = l 

See Figure [l] for a realization of ui n . 

We say that nodes v and w are connected at time n if there exists a path of 
distinct adjacent edges e%, ..,e m with w„(ej) = 1 for all i so that e\ connects v 
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Figure 1: Left: Input and output behavior at node v. The darkened line seg- 
ments indicate paths of sediment flow. The light line represents a potential flow 
path. Right: A small piece of a visualization of u n for some n. Only edges e 
with w„(e) = 1 are solid. 

to one of its children and e m connects w to one of its parents. Denote by C v>n 
the set of vertices which are connected to v at time n and define the backward 
(uphill) component of v = (x, y) at time n by 

C+ n :=C v>n n{( X ',y'):y'>y} 
Finally, let uj~ x = {e : uj n (e) = 1}. 

1.3 Regimes for 7] 

The parameter r\ plays an important role in the behavior of the model. For a 
fixed node v (at depth k) we have that for all n > 1 

lim F(d v {n) = L\d v (0) = R) = 

This indicates that when 77 is small the node v chooses a direction at time 
and has a high probability of sticking to this direction for most values of n > 1 . 
Since this is true for each node v, the evolution of {uj n } is somewhat simple. In 
the limit as 77 — > 0, each node picks a direction and stays with that direction for 
all time. That is, for each n > 1, and for each finite subset E C F? ev&n , 

iimP(c - 1 (i)ns- w - 1 (i)n^) = i (2) 

17 — >o 

and the dynamics has no effect on the configuration in any finite subset of 1? even . 
The configurations at any moment are the same as those in the discrete web 

m, m)- 

In the other direction, as r\ — > 00 each node "forgets" its history. That is, 
for each node v, the conditional probability given T n that it chooses left at time 
n + 1 is given in ([T]), and the limit of this quantity is 1/2. By symmetry, 
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P(d v (n + l) = L) = 1/2 = ¥(d v (n + 1) = R) 

and so 

lim [P(d v (n + 1) = L) - P(d v (n + 1) = L\T n )] = 0. 

rj — >oo 

Therefore the variables in any finite subset of {d v (n) : n > 0} converge in 
distribution to i.i.d. Bernoulli(l/2) variables. Furthermore, the variables in 
any finite subset of {d v {n) : n > 0,v <E l^ ven \ converge in distribution to 
i.i.d. Bernoulli (1/2) variables. Intuitively this holds because distinct nodes 
only interact with each other through their input and output loads, and both of 
these are eventually dominated by large ij. These statements indicate that when 
r\ is large, the dynamics of our erosion model are similar to those in a network 
in which each node flips a fair coin at each time n, independently from site to 
site and from time to time, to determine in which direction to send sediment. 
Thus the configurations {w n } resemble those taken from the dynamical discrete 

web am, 0). 

Given the relation both these extreme cases have to the variables {w„}, it is 
natural to view the present model (with < rj < oo) as an interpolation between 
the discrete web and the dynamical discrete web. Indeed, for each fixed n, the 
distribution of u n is the same as that in the case of rj = at time n = or that 
in the case of r\ — *■ oo at any time n. (In both cases, all directions are chosen by 
independent fair coin flips). 

As we shall see in section 3.2.2, the model with < r\ < oo can be likened 
to the case r/ — > oo in the following way. Each level k is associated to a measure 
9k (defined in eq. ([8|). Each node at level k samples (non-independently) from 
this measure a value p v . For any sequence n\,ri2, ■■■,n m of times and for any 
sequence x\, x m of elements from the set {L, R}, 

lim F(d v (n 1 +T)=x 1 ,...,d v {n m + T)=x m )=p^{l-p v ) NR 

T — >oo 

where Nl (Nr) is the number of i for which d{ = L. Because of this fact, we 
may view the model (for large time) as one in which each node fixes a Bernoulli 
parameter p v and flips a ^-biased coin independently each second n (but not 
independently from site to site) to determine the direction in which to move 
sediment. 

1.4 Outline of the Paper 

In section 2 we discuss the (relatively simple) behavior of nodes at depth 1. 
Since these nodes receive constant input load over time, we can use the well 
known model of Polya's Urn to analyze their output. In section 3 we discuss 
the more complicated behavior of nodes at depth at least 2. Here we make use 
of results of Pemantle [TJ for the time-dependent Polya Urn. We look more 
closely at properties of the input load, of the output load, and of the dynamics 
of these lower-depth nodes. 
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2 Top Level 



Since our top level nodes are equivalent to the model of Polya's urn, we recall 
basic facts of Polya's model Start with an urn containing Ro red balls and Bq 
black balls and draw one ball from the urn. Return this ball to the urn, along 
with another ball of the same color. After this round there are R\ red balls in 
the urn and B\ black balls in the urn, with either Rq = R\ or Bq = B%. Repeat 
this process infinitely many times, creating sequences {R n } n >o and {i? n }n>o so 
that for each n, 

R 

P(Rn+l — Rn = l\Rn,B n ) = — 

It is well known that the fraction — R R ™ B has an almost sure limit and 
that this limit is distributed as f3(Ro,B ) (see e.g. [8]). 

Let v be a node at depth 1. At the beginning of each second, v receives an 
amount of sediment equal to 1 and this input load amount does not change with 
time. The node sends this load either to the right or left, depending on the bias 
rule in . We are interested in the fraction of total load the node sends left 
(right) up to time n. To this end, define the load fraction 

LF v L (n) = , LF v «(n) = ^ , n > 1 

T v {n) J-v(n) 

Theorem 2.1. The quantities LF^{n) and LF^(n) have limits as n — ► oo. 
These limits are random: they are distributed as (3{rj 1 rj) . 

Proof. We will indicate the proof only for the case LF^(n). An easy calculation 
shows that Py(n) is a martingale w.r.t. {Tn} and, since it is bounded for all 
n, it has an almost sure limit. Solving for the limiting distribution is similar 
to solving for the related quantity in the standard Polya urn model. Sec, for 
instance, [5]. This gives 

rpL( \ i ^(") I g 

lim P v L (n) = lim ^ ( ") + j = lim *,(») + r.(n) 



ra— >oo T v (n) + 2?7 ri— >oo \ _|_ 



2r) 
T„(n) 



= lim = l im LF L( n) ( 3 ) 

n^oo l v [n) n—roo 

because r\ is constant w.r.t. n and T v (n) — > oo. □ 



Note that the limiting distribution in Theorem 2.1 is supported on [0,1] and 
has no atoms. This implies that with probability 1, the node v switches states 
(L,R) infinitely often and that neither of these states is transient. This is quite 
unlike the "sticking" associated to the dynamics in the 77 — > limit (refer to ([2])). 
The distribution from the above theorem for different values of r\ is pictured in 
Figure [2] For < 77 < 1 the limiting load fraction has a bimodal distribution, 
and for r\ > 1 the distribution is unimodal, symmetric about |. This means 
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0.25 0.5 0.75 1 0.25 0.5 0.75 1 0.25 0.5 0.75 1 



Figure 2: Asymptotic distributions for LF^(n): r\ — .5, 1,2 respectively. 

that when r\ is small, each node is likely to have a relatively strong preference 
for one direction and that when r\ is large, each node is likely to favor L and 
R somewhat equally. The case r\ = 1 gives a uniform distribution. Here v is 
equally likely to have a strong directional preference as it is not to. 

3 Lower Levels 

The simplicity of behavior at the top level comes from the fact that each node 
has an input load which is constant w.r.t. time. This is not true at lower levels. 
Each node has an input load whose magnitude is non-trivially time dependent. 
To make this more apparent, isolate an arbitrary node v at depth k. If at time 
t = n, v is not connected to either of its parents in u n , then its input load 
is 1 unit (coming only from rain). If, on the other hand, v is connected to 
at least one of its parents, then its input load will be strictly greater than 1 
unit. Therefore, the geometry of the connected components of u n determines 
the behavior of each node. This relationship is complex for at least two reasons. 
First, not only does the geometry of the network influence node behavior, the 
node behavior in turn determines the future geometry of the network. In this 
sense, our system generates its own randomness. Second, the method by which 
this randomness arises involves propagation. The geometry of nodes at depth 
k — I at time m affects the behavior of nodes at depth k at time n if and only 
if m = n — I. In other words, it takes I seconds for the output load from depth 
k — I to reach nodes at depth k. In spite of these complications, we set out to 
analyze these lower level nodes. 

The node v has an input load sequence I v = (I v (l),I v (2), ...), left out- 
put load sequence T„ — (Tjf (1), T^(2), ...), and output direction sequence 
d v = (d„(l), d v (2), ...). We are interested in analyzing the nature of this in- 
put sequence, the nature of the output sequence, and the relationship between 
the two. 

3.1 Input Load 

Figure [3] shows a histogram of input load values for all nodes at (a) depth 
5, (b) depth 7, and (c) depth 8 at t = 300s with r\ = 1 (the precise value 
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Figure 3: Load distribution for k — 5, 7, 8 respectively. 



of r\ does not matter, as a consequence of Theorem 3.1 1. The simulation was 
conducted with periodic boundary conditions, with 10 b nodes per row, and with 
10 rows. Therefore, the histogram for depth k at time n — 300s should closely 
approximate the probability mass function of the distribution of the input load 
for depth k at time n — 300s. One notices a few things. First, the support 
of the distribution at depth k is integers in the interval [1, \k{k + 1)]. Next, 
the mass function appears to decrease from load value 1 to a local minimum 
at k — 1, to increase for a bit to a local maximum, and then to decrease to 
the edge of its support. About 1/4 of nodes are at the heads of rills, while the 
fraction of rills starting short of the top increases with depth. The "bump" in 
the load distribution to the right of the value k—1 appears to travel to the right 
as depth increases. Looking at Figure [3] it is tempting to guess that the load 
distribution at a given level is a mixture of a distribution for loads that start at 
the top and one for loads that do not. Last, the different mass functions have 
several common values. For example, the probabilities for load values 1 to 4 are 
the same in each figure, and the probabilities for load values 1 to 6 are the same 
in the center and right figures. 

We present three structural theorems regarding the load distribution. The 
first gives basic information needed to make calculations, and the second gives us 
the value of the first moment of the distribution. The third discusses a limiting 
measure for the family of loads that do not originate at the top. Because of the 
simplicity of the first theorem, we state it without proof. 

Theorem 3.1. Basic properties of the load distribution. Let uq be a fixed time 
and let v be a node at depth k. 

a. All random variables d v (no) are i.i.d. with probability ^ of being L or R. 

b. The distribution of I v (no) is laterally translation invariant (i.e. along the 
x-axis) and is invariant in time for uq > k. 

c. The distribution of I v (jiq) is equal to the distribution of |C+ n | for any 
node w with depth equal to min(no,k). Therefore I v (no) takes values in 

^ n (7i + l) -| 
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Theorem 3.2. Let v be a node at depth k. The mean of the load distribution 
is 

Proof. We prove by induction on k. For k — 1, the statement is trivial, so 
consider fc > 1. Since the distribution of /„(n) is constant for n > fc, we assume 
n < fc. Let N Vj k-i be the number of nodes at level fc— 1 which send sediment to v 
at the end of time n — 1. This variable takes values in {0, 1, 2} with probabilities 
{1/4, 1/2, 1/4}, respectively. Call ui\ (1V2) the left (right) parent of v. 

2 

E(J„(n)) = ^E(/„(n)|JV„, fc _i = i)P(JV„, fc _i = i) 
i=l 



= 1/4 + 1/2 [1 + E(/ Wl (n - 1))] + 1/4 [1 + E(I Wl {n - 1) + I W2 {n - 1))] 



= 1 + E(J„, a (n - 1)) = 1 + (n - 1) = n 

where to go from the second line to the third line, we use the fact that the 
variables I Wl (n — 1) and I W2 (n — 1) have the same distribution (see b. under 



Theorem 3.1 1. □ 



Theorem |3.1| lets us use geometric properties of clusters of a static network 
(w no ) to study something which is dynamic: the load at time n at node v. That 
load may have come from a pathway that no longer even exists at time n. We 
further exploit this relationship, but to do this we must consider the concept 
of the dual web, defined in, for example, [7], and of whose definition we remind 
the reader. Consider the odd sublattice 



Kdd = {(x, y) G Z 2 : x + y odd and y < 1} 

For any node v* — (x* , y*) £ Z 2 dd we call the node (x* + 1, y* + 1) the right 
child of v* and we call the node (x* — l,y* + 1) the left child of v*. Similarly, 
we call the node (x* + l,y* — 1) the right parent of v* and we call the node 
(x* — l,y* — 1) the left parent of v* . We define the set E 2 ^ in the obvious 
way. The set of configurations {u> n : n > 0} induces a set of configurations 
{w* : n > 0} C {0,l} E ° dd by the following rule. If, in the configuration ui n , a 
node v = (x, y) is connected to its left child, we form a connection between the 
node v* = (x,y — 1) and its right child in the configuration by setting the 
image under w* of the edge in E 2 dd between v* and its right child to 1, and the 
image of the edge between v* and its left child to 0. If, on the other hand, v 
is connected to its right child in lo k then we set the image of the edge from v* 
to its left child under w* to 1 and the image of the edge from v* to its right 
child to 0. See Figure [3] and notice that we construct clusters in cj* so that no 
occupied edges in u/* cross any occupied edges in uj n . 
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Figure 4: Portion of an realization of the erosion network, along with its dual 
web. The solid lines indicate paths of sediment flow and the dotted lines show 
paths of the dual (courtesy of [B]). 

The upward paths in u>* now resemble the downward paths in u> n . That 
is, the upward path starting at the node v* is a simple symmetric random 
walk which is killed at depth 1. Random walks starting at different nodes are 
independent until they meet, at which point they coalesce into one random 
walk. (This is similar to the coalescing random walks picture of the discrete 
web, described in [TT]. 

There is an obvious physical interpretation for the paths in the dual web. For 
any two adjacent paths in the configuration u> n , there is a path in w* separating 
them. If the paths in u„ represent rills or drains, the paths in oj* represent the 
divides or ridges between them. Just as divides between rills do not cross rills, 
paths in w* do not cross paths in uj n . 

We now characterize the load distributions for our model. For any node 
v — (x, y) (with depth k), let v* L = (x — \,y) and let v* Ti — (x + 1, y). Consider 
the set of edges in the dual lattice contained in the paths emanating from the 
vertices v* R and v* L in Lu* t until cither (a) they meet at some vertex w* or (b) 
they reach a depth of 1. The set of nodes in 1? even in the interior of this set of 
edges is exactly the backward cluster of v in the configuration uj n . 

We now make some definitions so that we can work with this load distri- 
bution. Let {Xf 1 : i > 2} and {Xf 1 : i > 2} be independent sets of ran- 
dom variables (also independent of each other) which take the values 1 and -1 
each with probability \. For i > 2 let Y l = \{Xf - Xf) and for i > 1, let 
Wi = 1 + Y 2 + ... + Yi. Consider the stopping time 

r = min{ra : W n = 0} 
Up until the stopping time r, the random variable Wi represents the width 
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of the backward cluster of the node v in the real lattice (only valleys and not 
separating ridges), where we only consider nodes in this cluster whose depths 
are between k — i + 1 and k. Therefore the total number of nodes in this partial 
cluster should be 

L t :=W 1 + ... + W i 

Now we can make an equivalent definition of the distribution of the load 
I v (n) by saying that for each fixed n, it is the same as the distribution of the 
random variable 

Irfc(n) := £ m in(r,n,fc) (4) 

This variable is essentially a discrete integral of the symmetric random walk 
{Wi : i > 1}. 

Theorem 3.3. Let v be a node at depth k v and letw be a node at depth k w > k v . 
For any n > k v and for any I < k v we have 

V(I v (n) = l) = V(I w (n) = l) 

Therefore the limit 

lim I v (k v ) (5) 

k v — >oo 

exists in distribution. This limit is a.s. finite but has infinite mean. 
Proof. On the event r > k Vl i.e. the load originated from the top, 

Iv{ n ) = -^min(T,n,fc^) = ^kv > k v > I 

and 

Hence, we need only consider r < k v . 

P(/„(n) =l)= V(I v (n) =l,r< k v ) = P( J L min(r ,„ ) = 1,t < k v ) 

- P(imin(r,n,fe ro ) = l,T < k v ) = P(I w (n) = I) (6) 

The random variable Lfe u (n) is constant for n > k v , so 

F(I v (k v ) =l)= ¥(I v (k w ) =l) = ¥(I w (k w ) = I) 
where in the last equality we use (pj. Consequently, for any fixed I, the limit 

lim ¥(I v (k v ) = I) 

k v — >oo 
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exists. By the definition Q, a random variable with this limiting distribution 
is 



Lao '■— hm L min i T k k \ — L T 

k—>oo 

Since 

r<L r <±±ll 

the third statement of the theorem will follow if we show that r is a.s. finite and 
has infinite mean. But since the increments {Yi} of the random walk {Wi} have 
mean zero, the walk is recurrent. In addition, it is a standard result that the 
entrance time of the set {0} has infinite mean. This completes the proof. □ 



3.2 Dynamics 

Now we investigate some aspects of the effect of ?/ on the stability of config- 
urations over time. As noted, the dynamics creates an interpolation between 
the discrete web and the discrete dynamical web. The evolution of the system 
mirrors some aspects of rill erosion, one being that nodes through which a large 
amount of water passes at time uq have a non-trivial probability to channel a 
large amount of water at any time n\ > n . The degree to which this is true 
depends on the parameter n, as we will see. 



3.2.1 Load Correlation 

We start our analysis by inspecting simulation results. For any two positive 
integers M, N, let Vm,n be an enumeration of the MN nodes in the box [0, M — 
1] x [— N, —1] and define the load correlation coefficient at time n (for n > N) 
by 

v J2 v ev M , N K(N)I' v (n) 
KM,N{n) 



^/(E v ev M , N ^( N ) 2 )(j: ve v M , N lUn) 2 ) 

where I' v (n) = I v (n) — E(/„(n)). This quantity is only defined for n> N because 
two load vectors for a box of depth N are in some sense incomparable if they are 
taken at times no, ri\ with no < N < ni- For example, a node at depth n only 
has a maximum possible load of a t time n < N , whereas its maximum 

possible load is jV( -~^ +1 - ) for n> N. 

In Figure [5] we have graphed simulation results for a network 49 nodes deep 
and 49 nodes wide. The x-axis represents values of the parameter 77, as it 
varies from to 5. The y-axis represents values of a time averaged correlation 
coefficient, namely the quantity 

N' 

^ K M ,N{n) 



N' -N 
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1.2 




Figure 5: Averaged load correlation (vertical axis) versus n (horizontal axis). 
Data points are taken at ?y-intervals of 0.1 until rj = 2, and then at intervals of 
1. 



for M = N = 49 and N' — 200. Furthermore, we averaged this value over 6 
independent trials. This quantity is meant to approximate values of Km^n^) 
for M, N, n large. The time averaging seemed necessary because of fluctuations, 
most likely due to finite size conditions, in the quantity Km,n(ti). 

We can see in Figure [5] that the coefficient approaches 1 as n — > 0. This 
makes sense because, as remarked in section 1.3, the 77 — > limit of the dynam- 
ics (in any fixed box) is the same as the dynamics (or rather non-dynamics) of 
the discrete web. Therefore the load vector for this box should be similar (if not 
the same) at any two times. As r\ increases, the correlation coefficient decreases 
and appears to approach 0. Indeed, additional simulations give the following 
data: for r/ = 10, 100, 1000, 10000, the coefficients were .2391, .0896, .0189, and 
.0101,. From the discussion of the rj — > 00 limit given in section 1.3, the cor- 
relation coefficient should approach that computed from two load vectors from 
independent realizations of the discrete web. 

3.2.2 de Finetti Measures 

Whereas we can compare nodes at the top level to standard Polya urns, we can 
compare lower level nodes to time- dependent input [14] or random input |13j 
Polya urns. We start with an urn with Rq red balls and Bo black balls, as before, 
but we also have a time-dependent (or random) input sequence / = (Io,h, ■■■)■ 
At time t = n we draw a ball from the urn and we return it to the urn along 
with I n balls of the same color. Notice that this process with / = (1, 1, ...) is 
just the standard Polya urn. 

To analyze these lower level nodes, we will also make use of a fundamental 
result in the theory of exchangeable variables. 
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Definition 3.4. {0, l}-valued variables X\, X%, ... are exchangeable if for any 

X\, x m € {0, 1} and for any permutation a of m elements we have 

P(Xl = x x , -,X m = x m ) = P(X CT (i) = xi, ...,X a(jn ) = x m ) 

Theorem 3.5 (de Finetti). Let (O, J-",P) be a probability space and suppose that 
{X n } n >o are exchangeable {0, l}-valued random variables defined on fl. Then 
there exists a random variable F on £1 so that conditioned on F , the random 
variables X n are independent Bernoulli with parameter F. 



It is easy to verify that if depth(u) = 1, then the {0, 1}- valued var iabl es 
{D^(n)} n >o arc exchangeable. In our case, the variable F from Theorem 
actually 



3.5 



is 



p v := lim P v L (n). (7) 

n — >oo 

Thus if we know the asymptotic fraction of left choices for a node, then our 
node is just flipping independent coins each second with the same bias. 

At lower levels, the variables {D^(n)} are not exchangeable. However, they 
are asymptotically exchangeable. We use the definition of Kingman [T2"] . 

Definition 3.6. {0,l}-valued random variables Xi,X2,... are called asymp- 
totically exchangeable if there exists a sequence Yi,3^2, ... of exchangeable 
random variables so that for each x\, x m € {0, 1}, 

lim V(X 1+N = xi, ...,X m+N = x m ) 

./V— >oo 

= ¥(Y 1 =x 1 ,...,Y m = x m ) 



In the language of Theorem |3.5| let F be the random variable associated with 
the exchangeable variables {A„}. We call F the de Finetti measure for the 
sequence {Y n }. 

Let v be a node with depth k > 1. 

Theorem 3.7. The variables {P^ (n)} n >o form a bounded martingale sequence 
w.r.t. T n . Therefore they have an almost sure limit p v . 

Proof. Similar to the proof of (T3J Theorem 2.1] □ 
Rem ark 3.8. Using the same equations which produce ([3]), the limit in Theorem 



3.7 is the same as the limit of the variables {LF„ (n)} Ti 

For any number < p < 1, define the measure Q p on the set {0, 1} by 

Q p ({o}) = i-p,Qp({i})=p 

Let {«!, v r } be a finite set of vertices. For any vector of real numbers 
(pi, ...,p r ), each between and 1, define the product measure Qp on vectors 
in {0, l} r to be the product measure Y[l=i Qpi- 
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Theorem 3.9. For fixed v, the variables {D^(n) : n > 1} are asymptot- 
ically exchangeable with de Finetti measure equal to the distribution of p v . 
Furthermore, let v = (vi,...,v r ) be a vector of vertices and for each n, let 
D~(n) = (Z)J (n), D Vr {nj). If d±, d s are vectors in {0, l} r , with proba- 
bility one, 

s 

lim P(£>£(1 +T) = ds., ...,D$(s + T) = d s ) = E(]7 Q^dj)) 

T — >oo 

i=l 

where p = (p vi , ...,p Vr ). 

Proof. Similar to the proof of [T5] Theorem 2.2]. □ 

Because of lateral translation invariance, the de Finetti measure for v de- 
pends only on the depth k. In light of this, we define 

6k = de Finetti measure for row k (8) 

With this framework we will be able to study the switching rate of each node v 
once we have the following lemma. 

Lemma 3.10. For any node v, almost surely, 

1 - 

lim -r^(i)=p, (9) 

n— too n 

i=l 

Proof. Similar to the proof of [T3j Theorem 2.3]. □ 

We now define the switching function s v for n > 2 by s v (n) = D^(n)(l — 
D%(n - 1)) + D^(n - 1)(1 - D%(n)). Define the switching rate S v (n) to be the 
time average of s v , that is 



1 - 

S v (n) = } s v (i) 

n — 1 £ — ' 

i=2 



Theorem 3.11. The n — > oo limit of S v (n) exists a.s. 

lim S v (n) = 2p v (l —p v ) (10) 

n — too 

Proof. The proof is similar to the proof of [T3"l Theorem 2.3]. Let d n = D%(n) 
and p n = P^(n). A straightforward calculation gives 

E(d n+1 \T n ) = p n , E(p„ + i|.F n ) = p n (11) 

Now, 

^ n 1 n 

lim -S^ s v (i) — lim - ^(^(1 - di-i) + d»-i(l - dA) 

n — too fi * — * n — too Ji A — * 

i=2 i=2 
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^ n \ n 1 71 

= lim N + 2 lim — N didi-i = 2p v — 2 lim > cWi-i 

rj. — >r*i n ' ^ n. — >r>o T) ' ^ n — >oo J\ * 



by Lemma 3.10 We must show that the last limit above is a.s. equal to p\. 
Define M n = ^2™ =2 (di — p n )di_i. M n is a martingale with respect to T n 



E(M n+1 \J r n ) = E(y](djdj-i) - y^(p n+ idj_i) + d n+1 d n -pn+idnl^n) 



i=2 



i=2 



= } (djdi-i) - ^(pndi-x) + p n d n - p n d„ = M n 

i=2 (=2 

where we use both equations in ( |TT| ). Note also 

Y n 1 n 

lim - Pndi-i — p v lim -S^di-i 

n — >oo TL — ' n — >oo n * — • 



i=2 



i=2 



which is p 2 , , by Lemma 3. 10| Therefore it suffices to show that with probability 
one, 



Mr, 



0. 



(12) 



By summing the series 



- {M 1 + ... + M n = V — 
?! — * i 



i=2 



by parts, it can be shown that ( 12 1 will follow once we show that 

U.i M, 



E 

i=2 



converges. To this end, define M' n — Y^i=2 M ' +1 ^ M ' . We leave the reader to 
verify that M' n is a martingale. Using L 2 -orthogonality of martingale differences, 



i=2 



i=2 



" 1 

E ^2 



i=2 



Therefore is an L 2 bounded martingale and converges a.s. This completes 
the proof. 

□ 
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If p v € (0, 1) then neither of the choices L or R are transient for v. This 
prompts the question of whether or not the de Finetti measures 9k have atoms 
at or 1. For any fixed k, the answer is no. 

Theorem 3.12. For each k > 1, the measure 9k has no atoms. 

Proof. In |141 Theorem 4] it is shown that a time-dependent input Polya urn's 
de Finetti measure cannot have atoms if there is a C so that I v (n) < C for all n. 
For each realization of the dynamics and for each v, we have I v (n) < fc "( fc ^ +1 ) 
for all n. The result follows. □ 

Corollary 3.13. Each node v has a nonzero asymptotic switching rate. There- 
fore, for each v, the states L and R are recurrent. 

Proof. This is a direct consequence of Lemma [3T0| and Theorem [3T2| □ 

Corollary 3.14. With probability one, for each node v, the variable I v {n) takes 
each value in [1, fe "( fc ^+ 1 ) j j or infinitely many values of n. 

Proof. Let v\, ■■■,v m be the m — fct, ( fc ^ +1 ) — \ nodes above v which can send 
sediment to v and let d\, ...,d Vk £ {0, l} m . Let N > 1 and write Di(n) for the 
vector (£>£»,... ,D^Jn)). 



lim lim F(D§(T+jk v +l) = d[, D§(T+(j+l)k v ) = d k for some 1 < j < N) 

equals zero almost surely, by Theorem |3 . 9| and Theorem |3.12| But this probabil- 
ity dominates the probability of the event {Di [n + 1) = di , D~ (n + Vk + l) = 
d Vh for infinitely many n} c . The result follows. □ 

Here an interesting picture of our network emerges. On the one hand we may 
view the system as an infinite lattice (the lower half plane), where each node is a 
random input Polya urn. The output of the urns at depth k at time n becomes 
the input of the urns at depth k+1 at time n+1. On the other hand, as remarked 
in section 1.3, we may first sample (non-independently) values {p v : v € ^ ven \ 
from the de Finetti measures {9k : k > 1} to create an infinite array. As time n 
approaches infinity, the behavior of the system approaches the behavior of the 
same network in which each node v chooses to send its current load left with 
probability p v and right with probability 1 — p Vl independently at each second. 
Therefore this picture is of a network of two variables, a realization of values p v 
from the de Finetti measures, and realization of dynamics which coincides with 
the dynamics of a much simpler network. This second network is an obvious 
generalization of the Dynamical Discrete Web. 

Figure [6] shows histograms for the de Finetti measures 9k for k = 2, 5, 9 and 
for values of r\ = .5, 1, 2. One sees that the measures become more biased as k 
increases (for fixed rf). In other words, the mass of 9k is concentrated on domains 
closer to and 1 than is the mass of 9k-i- This would seem to imply that the 
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r\ = 1, row 2 77 = 1, row 5 77 = 1, row 9 




0.25 0.5 0.75 1 0.25 0.5 0.75 1 0.25 0.5 0.75 1 



77 = 2, row 2 rj = 2, row 5 77 = 2, row 9 




Figure 6: Defirietti measures for = 2,5,9 (from left to right) and 77 = .5, 1,2 
(from top to bottom). 



expected asymptotic switching rate of a node at level k (which is 2p v (l — p v )) 
must decrease with k. Similarly, if k is fixed and 77 decreases to 0, it seems that 
the expected switching rate should decrease. 

Figure [7] represents data given by simulations conducted with an erosion 
network with width 10 5 , depth 50, and 77 = .1, 1, or 10. The simulation ran for 
n = 1000 steps and at the end, switch rates for each node in the network were 
computed. In each row, each node's rate was averaged. Since two nodes v\ and 
V2 with the same depth k have independent behavior as long as they are at least 
a distance of 2k apart, the ergodic theorem gives that, as the network size ap- 
proaches infinity, the resulting average should resemble the expected switch rate 
for a row. The above results results were averaged by row over 3 independent 
trials. Finally, the data were plotted by row. Not only do the average switch 
rates appear to decrease as k increases, there appears to be a non-trivial (i.e. 
non-zero and 77 dependent) limit for the switch rate. This indicates that for at 
least some values of 77, the limit of the measures 9k (if it exists) is most likely 
not equal to |(<5 + <5i)- 
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Figure 7: Average switch rate (vertical axis) versus row (horizontal axis). The 
values of 77 are 10 (top curve), 1 (middle curve), .1 (bottom curve). 



3.2.3 Catastrophes 

Next, we study the following situation. Suppose a node v at a large depth k (for 
this section we assume the depth is at least 2) starts with a small input load 
and keeps a relatively small input load until a much later time. Then v's load 
changes dramatically. If this new load is sufficiently large, it could bring v's de 
Finetti measure much closer to h(So + Si). This analysis is from the point of 
view of the node v, whereas the analysis of the last half of the section will be 
from the point of view of the parent. 
Let 

A v (n) — , n > 1 

n 

Definition 3.15. For any n > 1, define the flood ratio F v (n) = ^ . For 
c > 1, we say that a flood of order c occurs at time n if F v (n) > c. 

Remark 3.16. Since I v (n),A v (n) g [1, ^(k(k + 1))], we have 

2 . . k(k+ 1) 

Proposition 3.17. For any v, A v (n) has a limit a.s. Therefore, 

liminf F v (n) < lim sup F v (n) 

n >oo n — 

Proof. We show the first statement by induction on depth(fc). Clearly this is 
true if depth(w) = 1. Otherwise, let wi be the left parent of v and assume that 
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for all nodes w' with depth equal to that of w\, A w i(n) has a lim it. Let N , R (n) 
be the number of i < n such that D R i (i) — 0. By Lemmas 



3.7 



and 



3.10 



T R (ri) T R (n)T w Jn) 

Km Wl = lim w \ Wl = (l-p Wl ) lim A Wl (n+l) 

n^oc n + L n^°o + L)l Wl (n) n^oo 

T L (n 

exists. The same argument shows that, if W2 is the right parent of v, linin^oo 
exists. Therefore, 

li m T ^ n ~ l ) = Hm n-l + T R (n-2) + T^(n-2) 

n — >oo fi n — >oo fi 

exists. 

For the second statement of the proposition, we use Lemma [3.14| to see that 

1 k(k + l) 

hmmf F v (n) = —j- -r < — —j- -r = hmsupF„(n) 

n->oo lim n _ 00 A„(ra) 21im n _ 00 A v (n) , woo 

□ 

Remark 3.18. The above proof shows also that linin^oo exists and equals 

(1 - p„) lim A v (n). 

Remark 3.19. A simple extension of the above proof using Theorem |3 . 1 2| shows 
that for any level k, the distribution of the time-average of the load for level k 
has no atoms. 

From the previous proposition we know that if v has an asymptotic average 
input load A v , then v will have infinitely many floods of order 2a"+e ^ or an y 
e > 0. This number can be quite large if A v is small. We study numerically the 
rates of these large floods as they relate to both depth(w) and 77. 

For a fixed node v, define the random measure 

n 
i=2 

Let {i>i}i>i be an enumeration of the nodes with the same depth as that of v. 
By the ergodic theorem, the average 

1 M 

1=1 

as M — > 00. See Figure [8] The graphs come from a simulation run for 10 5 
seconds on a network with a width of 10 4 nodes and a depth of 50 nodes. We 
graphed the density function for the measure j^g X)i=i ° [m^io 5 — Mu;,9ooo] in 
an attempt to approximate „_9 00 o E (/ i '",») f° r n l ar ge and for v with depth 5, 
20, and 50. The reason we subtracted fJ, Vi ,90oo is to decrease the effect of small 
times, during which the ratio F v (n) is likely to be an integer. 
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V = Jo, row 5 
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77 = jq- , row 20 
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77 = 1, row 5 



2 3 4 5 



3 4 5 



77 = 1, row 20 

It3.74 



4^ 



0.5 







0.5 



77 = 1, row 50 
T5.0 



012345 012345 012345 

77 = 10, row 5 77 = 10, row 20 rj = 10, row 50 



012345 012345 012345 



Figure 8: Histograms approximating the measure E(/i„.„) for n large. 



As 77 — > with a fixed row or as the row increases with fixed 77, each measure 
seems to concentrate its mass at and 1. In other words, the measure of any 
interval which does not include either of these two points appears to approach 0. 
In addition, Table[T]shows that as 77 — > with a fixed row or as the row increases 
with fixed 77, the expected fraction of time during which a large flood occurs 
(ratio above 5) increases. Since the time average of flood ratios approaches 1 as 
n — > 00, the above facts indicate a trend that as 77 — ^ or as the row increases, 
the time variance of measures increases, giving more possible variability of the 
flood ratios. Although this variance seems to increase, values near 1 (on the 
x-axis) show that the fraction of time flood ratios spend near 1 increases. In 
spirit, this is in accordance with previous results, as we explain. Figure [7] shows 
that as the row increases, the expected switching rate of a node decreases. It 
is reasonable to believe that the same conclusion holds if the depth is fixed but 
77 decreases. Therefore the network prefers to be more static in these circum- 
stances and we would expect a node to receive a relatively constant load, forcing 
flood ratios to be near 1. 

We now change focus to the parent. Make the following definition. 

Definition 3.20. For any n > 2, define the right catastrophe ratio 
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V 


Row 5 


Row 20 


Row 50 


.1 


.00035 


.00205 


.00308 


1 


.00136 


.01235 


.01789 


10 


.00024 


.02356 


.03711 



Table 1: Expected fraction of time during which a large flood (ratio at least 5) 
occurs. 



= "T^T whenever D^(n) = 

Here, A^{n) = yk?—^ > where N^(n) is the number of i < n such that D^(i) — 
0. For c > 1, we say that a right catastrophe of order c occurs at time 
n if C^(n) > c. Make similar definitions for left catastrophe ratio and left 
catastrophe. 



we 



Remark 3.21. From a similar argument to that used in Proposition 3.17 
have 

liminf C* (n) < limsupCf (n) 

We investigate the relationship between floods and catastrophes. Figure [9] 
shows the expected fraction of right catastrophes from the left parent which 
result in a flood of at least the same magnitude. The simulation was performed 
with a network of width 1000, depth 50, r\ equal to either .1, 1, or 10, and for 
a duration of 10 5 seconds. The calculation of fractions was only made between 
9000s and 10000s and we only consider catastrophes with ratio at least 1. It is 
clear from the figure that as the row increases, this expectation decreases. As 
rj decreases for a fixed row, the expectation also decreases. As n — > oo, 

^» = "> lim A ^ 

(n — l)l\^(n — 1) n-+oo 



by Remark 3.18 and Lemma 3.10| Therefore, we may use Remark 3.19 to show 



that almost surely for large n a right catastrophe of order c occurs for the node 
v at time n if v has a flood of order c at time n — 1. In other words, whenever 
a node receives a flood of order c at a large time, it has either a left or right 
catastrophe of the same order at the next second. 

Now we may interpret the probability that a node has a flood given that its 
left parent has a right catastrophe as the probability that a parent's catastrophe 
incites a catastrophe in the child. This would be a step of a possible catastro- 
phe cascade. The simulation results indicate that cascades become less present 
at lower levels (on average) but that they should never cease to exist. Two 
questions naturally arise. Given that a node has a right catastrophe of order c, 
how far does its catastrophe cascade travel? At each step of the cascade, the 
relevant (right or left) catastrophe ratio will generally increase. Indeed, a child 
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Figure 9: Expected fraction of right catastrophes from left parent which result 
in a flood of at least the same magnitude. The values of eta are 10 (top), 1 
(middle), and .1 (bottom). The data are plotted by row. 

node may even receive a catastrophe from both parents. How large does this 
ratio become in a typical cascade? 

4 Conclusion 

We have shown that the erosion model exhibits many properties of rill erosion. 
Each node chooses a random initial direction (right or left) in which to send 
sediment and further such choices become biased at a rate largely determined 
by the parameter eta. This is similar to the method by which rills are cut into 
a hillslope. As more water and sediment flows through a rill, a channel is cut 
deeper, giving reinforcement to the path, making it more likely to carry sediment 
in the future. Though the dynamics manifests itself through reinforcement, no 
fixed node can become fully biased (i.e. have a de Finetti measure equal to a 
sum of two delta masses). That is, since each node has a non-trivial asymptotic 
switching rate, sediment flow emerging from it will take both a left and right 
path a positive fraction of time. This rate of switching appears to decrease as 
we move further down the hill. 

There are a number of questions which deserve careful analysis. Do the 
measures Ok have a limit? If so, one would expect the limit to be the de Finetti 
measure associated with the "infinity process". To define this process, we start 
with a lattice of nodes which extends infinitely far in both positive and negative 
y directions. Since the behavior of a node v at time n in the present model 
depends only on the nodes in the n — 1 levels above it we may consider the 
input to the node v at time n in the infinity model to be a function of the 
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output of this finite number of ancestors. In the same way we have analyzed in 
this paper, it is possible to show that a dc Finctti measure for this process 
exists and that 

9^ = lim e n (n), (14) 

n — >oo 

where the term inside the limit is the measure given by 

6„{n)([a,b]) = ¥{P^(n) e [a,b]) for depth(u) = n 

Does the measure 6^ have atoms for some values of r/7 If so, is there a critical 
77* so that for < r\ < rj*, 8^ has atoms? If the limit of 9^ (assuming it 
exists) is not the same as #oo, does this limit have atoms and is there a critical 
77 associated with it? 
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